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The low energy electronic spectra of rotationally faulted graphene bilayers are studied using a 
long wavelength theory applicable to general commensurate fault angles. Lattice commensuration 
requires low energy electronic coherence across a fault and preempts massless Dirac behavior near 
the neutraUty point. Sublattice exchange symmetry distinguishes two famihes of commensurate 
faults that have distinct low energy spectra which can be interpreted as energy-renormalized forms 
of the spectra for the limiting Bernal and AA stacked structures. Sublattice-symmetric faults are 
generically fully gapped systems due to a pseudospin-orbit coupling appearing in their effective low 
energy Hamiltonians. 
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Coherent interlayer electronic motion in multilayer 
graphenes play a crucial role in their low energy prop- 
erties [T]. This physics is well understood for stacked 
structures with neighboring crystallographic axes rotated 
by multiples of 7r/3, including AB (Bernal), AA, ABC 
stackings and their related polymorphs Here the in- 
terlayer coupling scale typically exceeds 0.5 eV and pre- 
empts the massless Dirac physics of an isolated graphene 
sheet. Indeed experimental work on Bernal stacked bi- 
layers [HHl] demonstrates that their electronic proper- 
ties are radically different from those of a single layer 
[3 IB]. Yet, recent experimental work has revealed a 
family of multilayer graphenes that show only weak (if 
any) effects of their interlayer interaction. These include 
graphenes grown epitaxially on the SiC (0001) surface 
[5HTT|. mechanically exfoliated folded graphene bilayers 
[T^ and graphene flakes deposited on graphite [13] ■ A 
common structural attribute of these systems is the rota- 
tional misorientation of their neighboring layers at angles 
9 n7r/3. A continuum theoretic model has suggested 
that misorientation by an arbitrary fault angle induces 
a momentum mismatch between the tips of the Dirac 
cones in neighboring layers suppressing coherent inter- 
layer motion at low energy [T3j. In this interpretation, 
the Dirac points of neighboring layers remain quantum 
mechanically decoupled across a rotational fault |1H 1141 - 
[T5] accessing two dimensional physics in a family of three 
dimensional materials. 

This Letter presents a long wavelength theory of elec- 
tronic motion in graphene bilayers containing rotational 
faults at arbitrary commensurate angles. I find that 
the Dirac nodes of these structures are directly cou- 
pled across any commensurate rotational fault, produc- 
ing unexpectedly rich physics near their charge neutral- 
ity points. The theory generalizes previous approximate 
analyses (Hj by treating the lateral modulation of the 
interlayer coupling between rotated layers which is es- 
sential for understanding the low energy physics. Im- 
portantly, commensurate rotational faults occur in two 



distinct forms distinguished by their sublattice parity. 
Structures that are even under sublattice exchange (SE) 
are generically gapped (nonconducting) materials while 
those that break SE symmetry have two massive (curved) 
bands contacted at discrete Fermi points. Both these be- 
haviors derive from the the spectral properties of AA 
and Bernal stacked structures, and can be understood as 
energy-renormalized versions of these limiting cases. The 
gap in the faulted sublattice-symmetric states appears 
as a new feature specific to the faulted structures due 
to a pseudospin-dependence of the transmission ampli- 
tude across a twisted bilayer. These results provide the 
appropriate low energy Hamiltonian(s) for rotationally 
faulted bilayers superseding the massless Dirac model of 
an isolated sheet. 

The crystal structure of two dimensional graphene 
(Fig. 1) has a Bravais lattice spanned by two primitive 
translations ti = e~*'^/^ and ^2 = e"^^ with sublattice 
sites at = 0(l/\/3). We consider rotational stack- 

ing faults that fix overlapping A-sublattice sites at the 
origin and rotate one layer through angle with respect 
to the other, with translation vectors (^'1,^2) — e'^(ii,i2) 
and basis t'ai^b^ = ^^^ta{b)- A commensurate rotation 
occurs when T„j „ = mti -\-nt2 = m't[ +n't2 = T'^, ^, i.e. 
at discrete angles indexed by two integers m and n where 
e{m,n) = arg[(TOe-*''/6-hne"/6)/(ne-"/6_^^gi7r/6)]^ 

this notation AA stacking (all sites in neighboring layers 
eclipsed) has 9 — Q and Bernal stacking has 9 — 7r/3. 
Small angular deviations from the Bernal structure have 
indices m — 1 and large n. The \/l3 x \/l3 structures 
with 9 = 30° ± 2.204° structures observed by electron 
diffraction from epitaxial graphene on the Si (0001) face 
correspond to {m,n) — (1,3) and {m,n) = (2,5)Ji)J. 

Commensurate faults occur in two families determined 
by their sublattice exchange (SE) symmetry. With the 
A-sublattice sites at the origin, a commensuration is SE 
symmetric if B-sublattice sites are coincident at some 
other lattice position in the primitive cell. This occurs 
when TB+Tp^q = tb' +Tp, ,j, for integers (p, q) and {p' , q'), 
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FIG. 1: (left) Lattice structure of graphene with two sites in 
the primitive cell [A and B) and primitive translations ti and 
t2- (Right) Brillouin zones for the two layers in a rotational 
fault: the Brillouin zone corners labelled Km and K'^, are 
rotated by angle 9 to the points Km{0) and K'^iO) in the 
neighboring layer. 



requiring integer solutions io p — (m — n + Zniq)/[Zn). 
This occurs only when m — n is divisible by 3 and then 
the coincident B{B') sublattice site occurs at one of three 
possible threefold-symmetric WyckofF positions of the cell 
(e.g. r,„^„/3 in Fig. 2). The remaining Wyckoff positions 
are occupied by the A(A')-sublattice sites (the origin) 
and by overlapping hexagon centers (H,H'). When m — n 
is not divisible by 3 the only coincident site is the A-site 
at the origin, and the remaining two threefold symmet- 
ric Wyckoff positions are occupied by B-sublattice atoms 
of one layer aligned with the hexagon centers (H') of its 
neighbor. Rotational faults at angles 9 — 7t/3 — 9 form 
commensuration partners with primitive cells of equal ar- 
eas but opposite sublattice parities. Fig. 2 illustrates this 
situation for two partner commensurations at 30°— 8.213° 
((m, n) = (1, 2)) (left) and 30° + 8.213° ((m, n) = (1, 4)) 
(right). The limiting cases of Bernal (odd) and AA (even) 
stackings form the shortest period commensuration pair. 

Because of the rotation, the Brillouin zones of the 
two layers have different orientations (Fig. 1(b)) shift- 
ing their zone corners {Kj^, K'^) to rotated counterparts 
{Kjni9), K'„^{9)). The low energy electronic states of 
the decoupled layers have isotropic conical dispersions 
near each of these points with E{q) = ±/iWf^|q| where q 
is the crystal momentum measured relative to the cor- 
ner and vp is the Fermi velocity. These spectra are 
described by a pair of massless Dirac Hamiltonians for 
the K and K' points of the two layers [21j. Interlayer 
coupling is studied using a long wavelength theory that 
represents the low energy states as spatially modulated 
versions of the orthogonal zone corner Bloch states of 
the two layers, i.e. 4'(r) — J2a'^K,a{r)ua{'r). In the 
first star approximation, appropriate to the interlayer 
coupling problem, the basis states are the Bloch waves 
i^K,a = (l/\/3) s*^" retaining reciprocal lat- 

tice vectors that constrain the sum in 'ipK,a to the three 
equivalent corners of the Brillouin zone. The coupling 




FIG. 2: Geometry of the commensuration cells for bilayers 
faulted at 21.787° (left) and 38.213° (right). Red and blue dots 
denote the atom positions in the two layers. The white rhom- 
buses denote primitive commensuration cells with the same 
area for these two structures. The dashed yellow rhombus 
denotes a y/S x y/3 nonprimitive cell. The left hand structure 
is SE odd, with coincident atomic sites only on the A(A')- 
sublattice at the origin, the right hand structure is SE even 
with coincident sites on the A(A') and B(B') sublattices at 
threefold symmetric Wyckoff positions in the primitive cell 
and overlapping hexagon centers (H,H'). The density plot give 
the magnitude of the interlayer hopping potential discussed 
in the text. 



between layers is derived from an interaction functional 
of the form 



C/ = (l/2) / d2^T,(f)|*i(rO-*2(r)|' 



(1) 



which correlates the amplitudes and phases of the Bloch 
waves in neighboring layers. Here Tg is a periodic cou- 
pling potential with the translational symmetry of the 
commensuration cell. Using Eqn. (1) one finds that the 
interlayer interaction energy can be expressed in terms 
of the Fourier transforms of the slowly varying fields 
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X ^ {t{g)ulMu2Aq + K^m - Km' - G) + cc.) (2) 

e 

where N is the system size, Ao is the area of a graphene 
primitive cell and t{G) is the Fourier transform of the 
interlayer potential Ti{r) on the reciprocal lattice of the 
commensuration cell Q. 

The continuum theory of reference |14) is recovered 
from Eqn. 2 by retaining only its = terms, thus treat- 
ing the interlayer coupling as spatially uniform. In this 
approximation the states near the tip of the Dirac cone 
in one layer are coupled to three pairs of states at ener- 
gies ±W* = ±hvF\Km — Km{9)\ in its neighbor. At low 
energies, the effect of this coupling can be treated per- 
turbatively, preserving the Dirac nodes of two isotropic 
velocity-renormalized layer-decoupled Dirac Hamiltoni- 
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However, new physics arises from the G ^ contribu- 
tions in Eqn. 2 which mediate a direct couphng between 
the Dirac nodes of neighboring layers and prevent mass- 
less low energy behavior. To study it, note that the recip- 
rocal lattice of the bilayer is spanned by momenta with 
/owr integer indices^ = pGi+qG2+p'G[+q'G'2 [10]. Mo- 
mentum conserving couplings between K points in neigh- 
boring layers occur when Km — Km' {9) = Gip,q,p\q') 
with the angle 0(m, n) specified. This is an interlayer 
umklapp process where the spatial modulation of Ti^f) 
provides precisely the transverse momentum required to 
transport an electron between the Dirac nodes of neigh- 
boring layers. The momentum matching condition re- 
quires integer p solutions to p = (m — n)/Zn + qm/n^ 
and occurs only for supercommensurate structures with 
nonzero mod(m,3) = mod(n, 3). Importantly if this 
condition is not satisfied, momentum-conserving inter- 
layer couplings still occur, but instead through the anal- 
ogous intervalley umklapp process, i.e. K'^ — Km'{0) = 
Q{p,q,p' ,q'). These two possibilities are complementary 
and mutually exclusive: one or the other must occur if 
the rotational fault is commensurate. These two criteria 
distinguish SE-even and SE-odd structures, so that the 
SE-even structures require direct K — K{9) coupling and 
SE-odd structures K — K'{9) coupling. 

To understand the consequences of the interlayer inter- 
action one requires a theory for the Fourier coefficients 
t{Q) in Eqn. (2). These can be calculated from atomistic 
models, but their relevant properties are determined by 
symmetry. Note that the coupling function Ti{r) is a 
real periodic function with the translational symmetry 
of the commensuration cell. The structure function for 
the fi-th layer, n^(r) = J2me[i] J2a e'^^'."-(''-^" °) super- 
poses the six plane waves of the lowest star of recipro- 
cal lattice vectors G^^m producing a standing wave with 
maxima on atom sites and minima in hexagon centers. 
A useful model for the interlayer coupling potential is 
Tg{r) — Co exp[Cin{r)] where Ji{r) = rii +712 and Cq and 
Ci are constants; Ti is a superlattice-periodic function 
with maxima for coincident sites and with exponential 
suppression in regions that are out of interlayer registry. 
The grayscale plot in Fig. 2 show the spatial distribu- 
tion of Tg^r) where Ci is determined by matching the de- 
cay of the hopping amplitude between neighboring layer 
atoms as a function of small lateral offsets. This density 
plot shows that the interlayer amplitudes between ro- 
tated layers have coherent structures in the forms of five- 
fold rings (from overlapping misaligned hexagons) that 
are arranged to form two dimensional space-filling pat- 
terns. SE-odd structures are symmetric under threefold 
rotations while the SE-even structures retain a sixfold 
symmetry. The separable form Ti{r) = /i(r)/2(^ allows 
one to deduce a scaling rule for the Fourier coefficients: 

tiG) ~ iae-'>/^^/N,)J2^eii] hi.r^,)e-'^-^^ where the sum 
is over atomic sites in layer 1, a and b are constants. 



and Nc is the number of graphene cells (per layer) in the 
commensuration cell. For large Nc the prefactor decays 
as a power law of the cell size reflecting the fraction of 
atomic sites in good interlayer registry while the sum de- 
cays quickly as a function of Nc because of cancelling 
phases in its argument. 

The interlayer Hamiltonian can be expressed by a 3 x 3 
array of scattering amplitudes derived from the i(t/)'s 
giving the allowed transitions Km — > Km' {9). Three fold 
rotational symmetry requires that this matrix has the 
form 

Vps ^\ V2 Vo Vi ] 
\ Vi V2 Vq J 

where the pseudopotential coefficients Vi are are matrix 
elements of Ti. Completing the sum in Eqn. 2 projects 
this into the sublattice (pseudospin) basis and gives the 
2x2 interlayer transition matrices Hint seen by the Dirac 
fermions. The low energy Hamiltonian for an SE-even 
bilayer can be expressed as a 4 x 4 matrix (acting on the 
two sublattice and two layer degrees of freedom) 

and for the SE-odd bilayer 

where fT„ are Pauli matrices acting in the sublattice pseu- 
dospin basis of the n—th layer and vp is the renormalized 
Fermi velocity. The interlayer matrices Tij^j are 

^int = e-*^/2 ) ' " j 

"H^i shows that interlayer motion of an electron for SE- 
even faults requires a unitary transformation of its (A, B) 
sublattice amplitudes represented as an xy rotation of its 
pseudospin through angle ip. This angle is not defined ge- 
ometrically by the fault angle 9 but rather is determined 
by the relative magnitudes of the three pseudopotential 
matrix elements Vi. By contrast interlayer motion across 
a sublattice asymmetric fault involves only the ampli- 
tudes on its dominant (eclipsed) sublattice. The contin- 
uum model of [14] is recovered by setting Hint = 0. 

In either case, below an energy scale V the electronic 
spectra deviate from the massless Dirac form and inherit 
curvature from the interlayer coupling as shown in Fig. 
3. V 10 meV for commensurations at 61 = 30° ± 8.213° 
with Nc — 7 graphene cells per layer in their commen- 
suration cells. Nevertheless the forms of these spectra 
apply generally to any pair of commensuration partners. 
SE odd faults mix the degenerate Dirac bands gapping 
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FIG. 3: Low energy electronic spectra for SE-odd and SE- 
even faulted bilayers are illustrated using partner commen- 
surations at rotation angles Q = 21.787° (odd, left) and 
Q = 38.213° (even, right). These spectra are symmetric un- 
der rotations in momentum space. SE odd faults have massive 
bands that contact at Fermi points (left) and SE-even faults 
are gapped (right). The lower row gives the spectrum for a 
Bernal bilayer (left) and for an AA bilayer, which show related 
spectral properties. 



one pair on the scale V, leaving a second pair of mas- 
sive (curved) bands in contact at .E = 0. By contrast, 
SE even structures are fully gapped where the gap arises 
entirely from the pseudospin rotation in Eq. 5. Indeed 
for = these spectra consist of a pair of Dirac cones 
offset in energy by a bonding-antibonding splitting, and 
intersecting at £' = on a circle in reciprocal space. The 
pseudospin rotation lowers the symmetry of the bilayer 
and produces an avoided crossing of these states. For SE- 
even rotated bilayers the interlayer coupling describes a 
type of spin orbit coupling with the sublattice pseudospin 
index playing the role of the spin. 

These behaviors have precise analogs for the limiting 
cases of unfaulted Bernal and AA stacked layers. The 
Bernal bilayer (lower left Fig. 3) has exactly the struc- 
ture found for SE-odd faults but on an inflated energy 
scale ~ 0.5 eV, reflecting the full alignment of all atoms 
on a single sublattice. Similarly, for AA bilayers (lower 
right Fig. 3) interacting Dirac cones are displaced in 
energy but without a pseudospin rotation, so they inter- 
sect on rings at = 0. All intermediate commensu- 
rate site-centered rotational faults display either energy- 
renormalized Bernal-like or AA-like low energy spectra; 



the reduction of the energy scale is a measure of the loss 
of interlayer registry in the faulted bilayer. This corre- 
spondence can be deduced from the lattice symmetries 
of the density plots shown in Fig. 2. 

Since V < W* rotational faults open an energy "win- 
dow" V < E < W* in which the physics is well described 
by decoupled two-dimensional systems before their inter- 
layer coherence is apparent. The theory of small angular 
deviations from Bernal stacking [14] can be understood 
as a collapse of the energy scale V relative to W* . Yet, 
physics at the scale V is accessible to experimental probes 
and highly relevant to electrostatic gating and charge 
transport in structures derived from multilayer graphenes 
[5l [T2] particularly those with faults near 9 — n/6 found 
in epitaxial graphenes on SiC (0001) [IJL. Additionally, 
the many body physics of these systems depends cru- 
cially on the low energy structure of these spectra [H] 
and it can be expected to be quite different for different 
stacking sequences. Thus one can regard faulted mul- 
tilayer graphenes as presenting a family of new materi- 
als with properties that interpolate between single layer 
graphene and bulk graphite in an understandable and 
hopefully controllable way. 
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